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0.1  Technical  Objectives 

The  objective  of  this  project  is  to  develop  a  computer  aided  weld  design  tool  for  a  computer 
integrated  design  and  manufacturing  system. 

The  weld  design  tool  will  allow  for  the  determination  of  appropriate  welding  conditions  and 
if  needed  auxiliary  heating  for  a  specific  structural  design  such  that  welding  distortion  is  within 
tolerance.  The  database  will  be  integrated  with  Computer  Aided  Production  Engineering  (CAPE) 
software. 

The  project  contributes  to  a  collaborative  effort  by  Maglev  Inc  and  other  Universities  entitled 
’’Demonstration  of  Computer  Aided  Manufactuting  Techniques  for  the  Precision  Fabrication  of 
Large  Steel  Curved  Plate  Beam  Components  for  Shipbuilding  and  Other  Industries.” 

0.2  Technical  Approach 

The  proposed  program  is  organized  in  the  following  five  tasks: 

1.  Identification  of  structural  features  that  welding  distortion 

2.  Inverse  method  and  database  architecture 

3.  Generation  of  a  database 

4.  Integration  with  computer  aided  production  engineering  software 

5.  Modifications  and  refinement  using  production  tests 

0.3  Progress 

Progress  was  achieved  in  the  following  topics: 

•  Contributed  to  the  development  of  a  Lagrangial  thermo-mechanical  process  optimization 
approach.  The  methodology  was  demonstrated  in  the  optimization  of  the  thermal  tensioning 
process  for  the  minimization  of  welding  residual  stress  elimination  of  buckling  distortion. 

•  Developed  eigenvalue  and  large  deformation  sensitivity  analysis  capabilities  for  the  prediction 
of  buckling  distortion  and  bowing  distortion  in  large  welded  structures. 

•  Initiated  research  on  the  coupled  multi-scale  therm-mechanical  process  modeling  of  large 
structures.  The  approach  combines  adaptivity  and  domain  decomposition  for  both  temporal 
and  spatial  multi-scale  problems. 


Chapter  1 

Sensitivity  Analysis  and  Optimization 
of  Thermo-Elasto-Plastic  Processes  1 


1The  work  presented  in  this  chapter  has  been  submitted  for  review  to 
and  Engineering 


Computer  Methods  in  Applied  Mechanics 


20030418  094 
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1.1  Introduction 

1.2  Analytical  Formulation 

Finite  element  formulations  for  quasi-static  thermo-elasto-plastic  processes  in  Lagrangian  reference 
frames  have  been  widely  used  [1, 2,  3, 4,  5, 6, 7,  8, 9].  The  thermal  analysis  is  assumed  to  be  transient 
while  the  elasto-plastic  quasi-static.  Thermo-elasto-plastic  processes  are  typically  assumed  to  be 
weakly  coupled,  that  is,  the  temperature  profile  is  assumed  to  be  independent  of  stresses  and  strains. 
Thus  a  heat  transfer  analysis  is  performed  initially,  and  the  results  are  imported  for  the  mechanical 
analysis. 

1.2.1  Transient  Thermal  Analysis 

For  a  spatial  frame  x  fixed  to  the  body,  and  time  t,  the  governing  equation  for  transient  heat 
conduction  analysis  is  given  by  the  following: 

dT 

P^p~fc  =  ^  '  fcVT]  +  Q  in  the  entire  volume  V  of  the  material  (1.1) 

where  p  is  the  density  of  the  flowing  body,  Cp  is  the  specific  heat  capacity,  T  is  the  temperature,  k 
is  the  temperature-dependent  thermal  conductivity  matrix,  Q  is  the  internal  heat  generation  rate 
and  V  is  the  spatial  gradient  operator. 

The  initial  temperature  field  is  given  by 

T  =  T°  in  the  entire  volume  V  (1,2) 

where  T°  is  the  prescribed  initial  temperature.  The  following  boundary  conditions  are  applied  on 
the  surface: 

T  =  T  on  the  surface  AT,  with  prescribed  temperatures  T  (1.3) 

9  =  ?  on  the  surface  Aq,  with  prescribed  heat  fluxes  f  (1.4) 

Through  the  finite  element  formulation,  the  element  residual  vector  R  is  obtained  as  follows 

WJ 

+  ^,NTqwj  (1.5) 

Aq 

where  left  superscripts  n  -  1  and  n  represent  the  time  increments  n~H  and  nt,  T  is  element  nodal 
temperature  vector,  N  and  B  are  the  usual  matrices  which  interpolate  the  temperature  T  and 
temperature  gradient  VT  in  an  element;  J  and  j  are  the  volume  and  area  Jacobian  component 
corresponding  to  the  Gaussian  weighting  W  for  volume  and  w  for  surface  integration.  The  Global 
residual  vector  77  is  assembled  from  the  element  residual  vector  as  follows 

^(nr)=  E  RTT)  =  0  (1.6) 

Elements 

where  T  is  the  global  temperature  vector.  Equation  (1.6)  is  solved  in  an  incremental  iterative 
fashion. 


R(  ”T) 


=  E 

V 


BTkB  ”T  -  Nr<5  +  NtN pCl 


nT 


__  n— Irp  } 

-  n~H  j 
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1.2.2  Mechanical  Analysis 


V  •  S  +  b  =  0  in  V 


(1.7) 


where  S  is  tli©  second-order  stress  tensor ,  and  b  the  body  force  vector.  The  boundary  conditions 
are: 


u  =  u  on  surface  Au 

Sn  =  t  on  surface  A* 


(1.8) 

(1.9) 


where  u  are  the  prescribed  displacements  on  surface  Au,  t  are  the  prescribed  tractions  on  surface 
A*,  and  n  is  the  unit  outward  normal  to  the  surface  A*.  The  total  strain  is  the  Green’s  strain: 


E  =  \  {Vu  +  tVulT}  (1.10) 

Thanks  for  the  symmetry,  the  stress  tensor  S  and  strain  tensor  E  are  commonly  expressed  in  the 
vector  form  a  and  e  (usually  called  engineering  stress  and  strain)  for  the  computational  efficiency. 
Then  the  initial  conditions  are: 


(1.11) 

(1.12) 

(1.13) 


where  ep  is  the  plastic  strain  vector  and  eq  is  the  equivalent  plastic  strain. 

Assuming  small  deformation  thermo-elasto-plasticity,  the  total  strain  vector  e  is  decomposed 
into  the  elastic  strain  vector  ee,  ep  and  thermal  strain  vector  et: 


e  —  ee  +  ep  +  €t  (1*14) 

The  stress  strain  relationship  is: 

er  =  Cee  =  C  [e  —  €p  —  et ]  (1.15) 

where  C  is  the  temperature-dependent  material  stiffness  tensor. 

Through  the  finite  element  formulations,  the  element  residual  R  is  obtained  as  follows 

R(  «U)  =  £  [BT  V  -  NTbl  WJ-J2  NTt  wj  (1.16) 

v  A* 

where 

V  =  n-1«r  +  Act  (1.17) 

Differentiating  equation  (1.15)  using  incremental  expression  yields  following  equations 

Ao-  =  nC  [Ae  -  Aep  -  Aet]  +  AC  n_1ee  (1.18) 

The  elastic  predictor  aB  and  corresponding  elastic  strain  eBe  are  defined  as  follows 

=  n~lcr  +  nC  [Ae  -  Aet]  +  AC 
=  n  +  [Ae  —  Ae*] 


0\B 

«Be 


(1.19) 

(1.20) 
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Using  the  associative  J2  plasticity  [10],  the  yield  function  /  is: 

/  =  <rm  —  Oy  (1.21) 

where  crm  and  ay  are  the  Mises  stress  and  yield  stress.  Active  yielding  occurs  when  /  >  0.  In  case 
of  active  yielding,  the  evolution  of  Aeq  can  be  evaluated  by  the  radial  return  algorithm  [1], 


Aep  =  Aeqa 

(1.22) 

where 

a  = 

lLm 

(1.23) 

m  — 

1 

_  SB 
<*Bm 

(1.24) 

L  - 

diag(  1  1  1  2  2  2  ) 

(1.25) 

where  <TBm  and  s b  are  the  Mises  stress  and  the  deviatoric  stress  of  the  elastic  predictor  erg. 


1.3  Sensitivity  Analysis  from  the  Method  of  Direct  Differentiation 

1.3.1  Thermal  Sensitivity 

Differentiating  Equation  (1.6)  with  respect  to  each  of  the  design  variable  fa  yields: 


tPT 

d(j)% 


dn 


dnT 


J  constant  <f> ■ 


i 


vn 


(1.26) 


where  V  stands  for  the  contribution  due  to  all  variables  except  for  T.  Because  the  global  stiffness 
has  already  been  assembled  to  evaluate  the  global  temperature  and  every  global  expression  is  the 
assembly  of  elemental  expression,  the  only  term  that  needs  to  be  evaluated  for  thermal  sensitivity 
is  Differentiation  of  Equation  (1.5)  yields: 


VR 

V<(>i 


e{ 


dB 

9<f>i 


dB 


kB  nT  +  Bxk^-  nT  +  Br^B  nT  -  N' 


d<f>, 


9<f>i 


-dQ 

9<t>i 


r  T  -1 
d<j>i  H  -  n~H  +  P  p  H  -  n-it 


+ 


BTkB  ”T  -  N tQ  +  NrNpC3 


lp  nt  -  n~lt 


w 


WJ 

dj 


+E 


!SiT-~wj  +  N  Tqw 
o<Pi 


di 


3 


(1.27) 


1.3.2  Mechanical  Sensitivity 

In  this  section,  it  is  assumed  that  temperatures  and  their  sensitivities  are  available  (Equation 
(1.27)).  Further,  mechanical  forces  are  not  considered,  that  is,  temperature  change  is  the  only 
loading  in  this  sensitivity  formulation.  The  left  superscript  n  is  dropped  for  simplicity. 
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Primary  Sensitivity 

Differentiating  the  global  residual  equation  with  respect  to  each  design  variable  fa  yields 


my1 

dJU  „„ 


constant  (j>i 


(1.28) 


where  D  stands  for  the  variance  due  to  all  variables  but  U.  Similar  to  the  thermal  sensitivity 
formulation,  the  only  term  that  needs  to  be  evaluated  for  displacement  sensitivity  is  which  is; 


DR  ^  dBT  T  ~  D<r  «  dJ ' 

=  Vhwj+B  ,rww]  (U9) 

In  case  of  non-active  yielding: 

Va  Va  b 

vfa  =  d&~  (L30) 

In  case  of  active  yielding,  a  can  be  expressed  as  follows  from  the  radial  return  algorithm 

cr  =  crh  +  aym  (1.31) 

where  cr^  is  the  hydrostatic  stress  of  a.  Then  ^  is  evaluated  as  follows 


(1.30) 


i Vcr  —  Vcr},,  (  Dm  i  _ Doy 

Vfo~Vh+aYVh+mVfc 


(1.32) 


Secondary  Sensitivities 

The  secondary  sensitivities  can  be  evaluated  only  after  the  displacement  sensitivity  is  evaluated. 

T _ _ _  r _ _ ,  * _ _ *  l  j*  i 


In  case  of  non-active  yielding,  erg  and  ege  are  equal  to  a  and  ee. 

da  _  das  _  dasd U  (  Derg  „  ooX 

dfi  ~  d^"  aUd^+  D&  ^L33^ 

_  d,€.Be  ^eBe  dXJ  'D(,£e 

d<j>i  ~  d<f>i  ~  dv 

deQ  dn~1e„ 

dfa  ~  ~dfa~  ^L35) 

In  case  of  active  yielding,  the  secondary  sensitivities  are  evaluated  as  follows 

da  _  da  dU  Va 

dfa  ~  aU  d$i  +  Vfa  ( 1,36 ) 

dee  _  dee  dXJ  Dee 

dfa  ~  W  dfa  +  Vfa  ^,37' 

de„  de„  dU  Ve„ 

#i  “  aud^  +  D^  ^L38) 

Following  relation  can  be  obtained  from  Equation  (1.20),  Equation  (1.22)  and  Equation  (1.24): 

3 

€e  =  €Be  -  Aep  =  €Be  ~  Aeg-Lm  (1.39) 

Using  Equation  (1.31)  and  Equation  (1.39),  the  unknowns  in  Equation  (1.36),  Equation  (1.37) 
and  Equation  (1.38)  can  be  evaluated. 


In  case  of  active  yielding,  the  secondary  sensitivities  are  evaluated  as  follows 


da 

da 

dU 

f 

Va 

d(j)i 

au 

d<f>i 

T 

D& 

dee 

dee 

dU 

+ 

Dee 

d(bi 

au 

d<t>i 

deq 

deq 

dU 

4. 

Deg 

d<pi 

au 

d(pi 

T 

D* 

£e  =  €Be  ~  A £p  =  €ge  ~  Aeg-Lm 


8 


1.4  Numerical  Implementation 

Welding  causes  the  material  to  have  a  permanent  distortion  and  residual  stress  [11,  12,  8,  13]. 
The  transient  thermal  tensioning  also  known  as  side  heating  technique  can  be  used  to  control  the 
welding  distortion  and  residual  stress  without  modifying  design  specification  [14,  15].  There  are 
many  design  variables  that  characterize  the  side  heaters  such  as  the  heat  source,  position  and  shape 
of  side  heaters.  In  this  section,  the  heat  source  and  positions  of  side  heaters  are  optimized  with 
other  variables  fixed  for  the  minimum  residual  stress  using  the  sensitivity  equations  developed  in 
the  previous  section  in  a  3D  Lagrangian  reference  frame.  No  constraints  are  considered  in  this 
example. 

1.4.1  Welding  Conditions 


Welding 

direction 


Side  heaters : 

*  Moving  along  with  the  torches 


Figure  1.1:  Configuration  of  welding  and  side  heating  setup. 


The  schematic  welding  configuration  in  this  simulation  is  shown  in  the  Fig  1.1.  Side  heaters 
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are  followed  by  two  welding  torches.  Convection  boundary  conditions  are  assigned  for  all  free 


Figure  1.2:  Conductivity  (k),  specific  heat  ( Cp ),  and  air  convection  ( h )  for  A36. 


surfaces.  The  internal  heat  generation  rate  by  the  welding  torch,  modeled  with  a  ’’double  ellipsoid” 
heat  source  model  [16],  is  given  as, 


Q  = 


&V3Qwr)wf 

aboTsfir 


[W/mm3]  (1-40) 


where  Qw  (2680.35  W/mmB)  is  the  welding  heat  input;  rjw  (1.0)  is  the  welding  efficiency,  x,  y,  and 
^  are  the  local  coordinates  of  the  double  ellipsoid  model  aligned  with  the  weld  fillet;  a  (5\/2  mm)  is 
the  weld  width;  b  (5y/2  mm)  is  the  weld  penetration;  c  is  the  weld  ellipsoid  length;  v  (6.35  mm/s) 
is  the  torch  travel  speed.  The  numbers  in  the  parentheses  are  the  values  which  are  used  for  this 
implementation.  Goldak  et  al.  [16]  used  c  =  a  and  /  =  0.6  before  the  torch  passes  the  analysis 
region,  and  c  —  4a  and  /  =  1.4  after  the  torch  pases  the  analysis  region.  However  in  this  paper, 
c  =  4a  and  /  =  1.0  are  used  instead  to  improve  the  convergence  in  the  simulation.  In  fact,  these 
factors  have  a  measurable  effect  on  the  temperature  field.  However  they  have  negligible  effect  on 
the  residual  stress.  The  side  heat  source  is  applied  on  the  top  surface  of  the  plate  as  shown  in  Fig 
1.1  and  is  defined  as  follows 

q(x,z) 

Mx 

Mz 

where  x  and  z  are  the  local  coordinates  from  the  center  of  the  side  heating;  Qs(W/mm2)  is  the 
side  Heating  input,  t?s  (1.0)  is  the  side  heating  efficiency;  Bs  (6”)  and  Ts(l”)  are  the  band  width 
and  length  of  the  side  heating,  Sxi  (0.2),  Sx2  (0.2),  Szl  (0.2)  and  Sz2  (0.2)  are  used  to  control 
the  gradient  of  heat  flux  in  the  side  heater  edges.  This  side  heater  shape  is  shown  in  Figure  1.4. 


-  d-41) 

=  {tanh  (Sx2[a:  +  fa  +  -Bs/2j)  -  tanh  (Sxl  [x  +  fo-  BJ 2]) 

+  tanh  (5*i [®  -<f>2  +  Baf 2])  -  tanh  (Sx2[a:  -<h~  BtJ 2])}  /2  (1.42) 

=  {tanh  (52l[z  -  Ls/2\)  -  tanh  (Sz2[z  +  Ls/2])}  /2  (1.43) 
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T  (*C) 


Figure  1.3:  Elastic  modulus  ( E ),  Poission's  ratio  (u),  thermal  expansion  coefficient  (a)  and  yield 
strength  (ay)  for  A36. 


The  numbers  in  the  parentheses  are  the  values  which  are  used  in  this  simulation.  The  material 
properties  of  A36  steel  used  in  this  simulation  are  shown  in  Figure  1.2  for  the  thermal  analysis  and 
in  Figure  1.3  for  the  mechanical  analysis.  The  isotropic  hardening  coefficient  is  assumed  to  be  8000 
[MPa]  at  any  temperature. 


i 

o.s 

0.4 

0.2 

q 


i 

0.8 

5N0.6 

0.4 

0.2 

Q 


150 


Figure  1.4:  Side  heater  shape  parameters  Mx  and  Mz  (  see  Equation  (1.42)  and  Equation  (1.43)) 

A  finite  element  model  is  developed  as  shown  Figure  1.5.  The  dimensions  are  12”  x  12”  x  1/8” 
for  base  plate  and  12”  x  2”  x  1/8”  for  the  stiffener.  This  model  has  13864  nodes  and  2352  20-noded 
brick  elements.  Since  high  temperature  gradients  are  prevalent  at  the  welding  region,  the  mesh 
is  finer  along  the  welding  torch  path  and  coarser  away  from  it.  The  boundary  conditions  for  the 
mechanical  analysis  are  shown  in  Table  1.1. 
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1.4.2  Optimization 

Since  the  residual  longitudinal  compressive  stress  away  from  the  weld  zone  can  be  used  as  a  criterion 
of  welding  induced  buckling  [17],  the  optimization  problem  is  expressed  as  follows 


min  F  =  min  ^2(lecrezz)2  (1.44) 

e 

(*)  max  (1.45) 

where  a\z  is  the  centroid  longitudinal  residual  stress  at  element  e  in  the  objective  region  shown  in 
Figure  1.1,  and  le  is  the  x-direction  length  of  the  element.  The  gradient  of  the  objective  function 
F  is  obtained  as  follows. 


*T  -  2  )  oz  — 


(1.46) 


The  Design  variables  are  the  side  heat  source  Qs(=  <pi),  transverse  position  of  the  side  heater  <f>2 
and  the  distance  between  the  side  heater  and  the  first  welding  torch  as  shown  in  Equation  (1.41) 
and  Figure  1.1. 

The  optimization  loop  is  implemented  using  the  BFGS  line  search  method  provided  in  the  DOT 
package  [18]. 

Thermal  and  mechanical  analyses  and  their  sensitivity  analyses  are  performed  in  an  in-house 
SMP  FORTRAN  90  code. 


1.4.3  Results 

The  results  of  the  numerical  optimization  are  summarized  in  Table  1.2.  The  total  analysis  time 
for  each  side  heating  configuration  is  set  up  to  3000  seconds  for  both  the  thermal  and  mechanical 
analyses  because  temperature  distribution  becomes  uniform  and  the  residual  stress  distribution 
shows  no  more  change  after  that  time.  Since  a  transient  analysis  has  many  incremental  results, 
the  figures  shown  in  this  paper  are  selected  at  only  one  increment.  The  result  plots  of  temperature 
analysis  are  chosen  when  all  the  heat  sources  appear  in  the  analysis  model  and  those  of  mechanical 
analysis  at  the  final  increment. 

The  temperature  profile  at  the  initial  design  point  is  shown  Figure  1.6.  The  sensitivity  of 
the  temperature  field  with  respect  to  the  side  heat  source  at  the  initial  design  point  is  shown  in 
Figure  1.7.  Increase  in  side  heat  source  will  result  in  temperature  increase  in  the  side  panel.  The 
longitudinal  residual  stress  profile  at  the  initial  design  point  is  shown  in  Figure  1.8.  The  stress  is 
tensile  near  the  welding  torch  path  and  compressive  away  from  it.  Figure  1.9  shows  the  longitudinal 
residual  sensitivities  with  respect  to  the  side  heat  source.  The  stress  in  the  objective  region  (see 
Figure  1.1)  is  compressive  as  shown  in  Figure  1.8  so  that  positive  sensitivity  is  the  desirable  direction 
for  the  minimum  residual  stress  in  that  region. 

The  longitudinal  residual  stress  field  with  the  optimum  design  variables  is  shown  in  Figure 
1.10.  Compared  with  Figure  1.8,  the  residual  stress  reduction  is  observed  not  only  in  the  objective 
region  but  also  over  the  outside  panel.  Figure  1.11  shows  the  longitudinal  residual  stresses  along 
the  ’’Center  Line”  (see  Figure  1.1)  for  three  cases.  The  residual  stress  in  the  objective  region  is 
successfully  reduced  for  the  optimum  side  heater.  The  vertical  dotted-line  in  the  left  side  from  the 
axis  line  indicates  where  the  objective  region  start. 

The  variation  of  the  objective  function  defined  in  Equation  (1.44)  during  this  optimization  is 
shown  in  Figure  1.12.  A  total  of  28  function  calls  and  6  gradient  calls  were  made  during  entire 
optimization. 


Figure  1.9:  Sensitivity  of  longitudinal  residual  stress  with  respect  to  side  heat  source  at  the 
initial  design  point  [ MPa/W ] 
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2.1  Introduction 

Due  to  the  advantages  in  design  flexibility,  cost  savings,  reduced  overall  weight  and  enhanced 
structural  performance,  welding  has  been  applied  increasingly  in  comparison  with  other  mechanical 
joining  processes.  However,  several  types  of  distortion  can  be  induced,  as  discussed  in  detail  by 
Masubuchi  [19]  (Figure  2.1),  with  the  application  of  welding  process. 


Transverse 

Shrinkage 


Angular  Rotational 

Change  Distortion 


Longitudinal 

Shrinkage 


Buckling 

Distortion 


Longitudinal 

Bending 


Figure  2.1:  Types  of  welding  distortion  [19]. 

Among  the  types  of  distortions,  the  out-of-plane  deformations,  such  as  bending,  buckling 
and  angular  deformations,  significantly  influence  the  required  dimensional  precision  in  structural 
components.  In  order  to  acquire  reduction  in  overall  weight  and  achieve  more  controllable 
manufacturing,  relatively  thin  components  made  of  higher  strength  steels  are  preferred  when 
fabricating  large  structures.  However,  welding  in  thinner  components  introduces  a  significant 
drawback,  namely  out-of-plane  distortions.  These  types  of  distortions  cause  loss  of  dimensional 
control,  structural  integrity,  fit-up  between  panels,  and  thus  increase  fabrication  costs. 
Consequently,  analysis  on  the  out-of-plane  distortion  effects  is  important  in  many  industries  such 
as  shipbuilding,  railroad,  aerospace,  and  mass  rapid  transportation  systems. 

Finite  element  techniques  have  been  used  in  the  prediction  of  welding  residual  stress  and 
distortion  for  more  than  two  decades.  Due  to  the  nature  of  the  process,  additional  complexities  are 
involved  in  the  FEA  of  welding  compared  to  traditional  mechanics,  such  as  temperature  and  history 
dependent  material  properties;  high  gradients  of  temperature,  stress  and  strain  fields  with  respect  to 
both  time  and  spatial  coordinates;  large  deformations  in  thin  structures  and  phase  transformation 
and  creep  phenomena.  Welding-induced  buckling  of  thin-walled  structures  has  been  investigated 
in  greater  detail  in  [20,  21,  22].  Ueda  et.  al.  [20,  23]  presented  a  methodology  to  determine 
the  buckling  behavior  of  plates  by  large  deformation  elastic  FEA  and  employing  inherent  strain 
distributions.  Michaleris  et.  al.  [21,  24]  developed  a  predictive  buckling  analysis  technique  for  thin 
section  panels,  combining  decoupled  weld  process  simulations  and  eigenvalue  buckling  analyses. 
Tsai  et.  al.  [22]  studied  the  distortion  mechanisms  and  the  effect  of  welding  sequence  on  panel 
distortion.  No  published  work  is  available  for  the  computation  of  the  combined  effects  of  both 
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bending  and  buckling  distortion. 

Miehaleris  et  al.  [25],  demonstrated  that  3D  finite  element  models  of  the  welding  process  are 
needed  to  accurately  compute  angular  distortion.  Furthermore,  3D  weld  process  finite  element 
models  can  easily  account  for  the  increased  stiifness  of  plate  curvature  and  the  compliance  of 
fixturing  restrains.  However,  3D  finite  element  simulations  of  welding  large  structures  require 
very  large  models  where  heat  balance  and  equilibrium  is  iteratively  computed  for  several  thousand 
increments  [26],  Such  simulations  would  require  prohibitively  costly  numerical  computations  using 
currently  available  software.  The  development  of  an  efficient  computational  approach  is  needed  for 
performing  direct  3D  welding  simulations. 

In  this  work  the  decoupled  2-D  and  3-D  finite  element  analysis  technique  by  Miehaleris  et. 
al.  [24,  21]  is  applied  to  evaluate  welding-induced  longitudinal  bending  and  buckling  in  various 
stages  of  the  fabrication  of  the  maglev  guideway  beam.  Longitudinal  residual  stress  effects  are 
considered  only,  and  the  angular  distortion  which  usually  has  small  magnitude  [17]  is  neglected.  It 
is  demonstrated  that  1)  eigenvalue  analyses  can  only  compute  the  critical  buckling  stress  that 
induces  the  buckling  distortion  as  well  as  the  various  buckling  modes,  2)  small  deformation 
analyses  only  compute  the  bending  effect  on  the  structure  and  the  deformation  magnitude,  and 
3)  large  deformation  analyses  can  compute  the  combined  effects  of  bending  and  buckling  and  the 
deformation  magnitude. 

The  Maglev  system  is  a  high-speed  magnetically  levitated  ground  transportation  system  that  is 
designed  to  operate  at  speeds  in  excess  of  310  mph  at  a  high  level  of  ride  comfort.  The  system  has 
no  wheels  or  moving  parts  and  is  levitated  and  propelled  by  a  long  stator  linear  motor  embedded 
in  the  guideway.  Both  propulsion  and  levitation  are  provided  by  electro-magnetic  waves  resulting 
in  a  system  that  moves  contact  free,  with  no  wear  and  tear  on  the  system  and  no  friction  to  impede 
its  efficiency.  In  order  to  achieve  high  ride  comfort  at  speeds  in  excess  of  310  mph,  guideway 
structure  must  be  manufactured  within  very  small  tolerances.  The  47-mile  proposed  Pennsylvania 
alignment  consists  of  over  2000  guideway  beams,  each  measuring  203  feet  long,  weighing  135  tons 
with  compound  curves  built-in  and  having  to  be  manufactured  within  millimeters  of  tolerance. 
Precision  fabrication  technology  needs  to  be  developed  for  the  production  of  the  guideway  beam 
within  specifications. 

Figure  2,2  shows  a  section  of  the  maglev  guideway  beam,  which  is  the  main  component  of  the 
magnetic  levitation  transportation  system.  The  guideway  beam  is  double  span  and  is  supported 
by  piers  with  varying  distances  between  them  depending  on  the  beam  type  and  curvature.  The 
section  of  a  beam  known  as  the  Type  1  guideway  beam  is  analyzed  in  this  work.  As  shown  in 
Figure  2.2,  the  guideway  beam  is  a  trapezoidal  box  beam  structure  with  25  mm  thick  stiffeners 
located  at  fixed  intervals.  The  main  components  are  the  top  flange,  also  known  as  the  ’deck  plate’, 
which  is  18  mm  thick,  the  side  web  plates  which  are  12  mm  thick  and  the  bottom  flange  which  is 
30  mm  thick.  The  top  flange,  side  web  plates  and  bottom  flange  are  welded  longitudinally  using 
fillet  welds.  The  stiffeners  are  welded  onto  the  top  flange  using  double  fillet  welds. 

Dimensional  accuracy  is  especially  important  in  maglev  guideway  beam  technology  to 
achieve  high  level  of  ride  comfort.  Therefore,  the  out-of-plane  distortions  resulting  from 
welding,  such  as  bending,  buckling  and  angular  deformations,  need  to  be  kept  at  a  minimum. 
Furthermore,  if  buckling  distortion  occurs,  the  structural  instability  resulting  from  the  out-of¬ 
plane  distortion  violates  the  tolerance  requirements.  This  results  in  either  reduced  quality  in 
passenger  transportation,  or  even  re-manufacturing  of  the  structure  which  detrimentally  increases 
manufacturing  cost. 
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Figure  2,2:  The  components  of  the  maglev  guideway  beam 


2.2  Analysis  Approaches 

Following  the  work  of  Michaleris  et.  al,  [24,  21],  the  response  of  the  maglev  guideway  beam  is 
evaluated  in  two  steps  by  combining  three-dimensional  structural  analyses  with  two-dimensional 
welding  simulations  in  a  decoupled  approach.  The  longitudinal  residual  stresses  only  are  considered, 
while  the  angular  distortion  is  neglected, 

2- D  Thermo-mechanical  Weld  Simulation  : 

A  two  dimensional  thermo-elasto-plastic  analysis  is  performed  to  determine  the  residual  stress, 
and  plastic  strain  fields  during  the  welding  process.  The  longitudinal  residual  stresses  are  caused 
by  the  negative  plastic  strains  resulting  from  the  welding  thermal  cycle. 

3- D  Structural  Deformation  Simulation  and  Eigenvalue  Analysis  : 

The  longitudinal  bending  (bowing)  and  buckling  distortion  effects  are  consequently  determined 
by  applying  the,  mostly  uniform  and  compressive,  longitudinal  plastic  strain  field  of  the  2-D  weld 
model  on  the  3-D  structural  model  as  equivalent  load. 

A  constant,  negative  thermal  load  is  applied  at  the  weld  region  to  introduce  the  effects  of 
welding  into  the  3-D  structure.  Thermal  loading  is  used  rather  than  mapping  the  plastic  strain 
field,  which  would  require  a  complex  mapping  procedure. 

The  deformation  analyses  are  performed  with  applications  of  linear  (small  deformation)  and 
non-linear  (large  deformation)  formulations,  respectively.  With  the  linear  structural  deformation 
analysis  only,  the  bowing  distortion  magnitude  induced  by  applying  the  welding  can  be  evaluated. 
Meanwhile,  the  combined  bowing  and  buckling  distortion  magnitude  can  be  shown  in  the  non-linear 
results. 
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An  eigenvalue  analysis  is  performed  to  determine  the  critical  residual  stresses.  The  various 
buckling  distortion  modes  can  also  be  attained. 

2.2.1  Welding  Simulation 

The  welding  simulation  involves  a  thermal  and  a  mechanical  analysis.  The  effect  of  mechanical 
response  is  assumed  to  be  negligible  on  the  thermal  behavior,  thus  the  temperature  field  is  solved 
independently  from  the  mechanical  solution.  To  determine  the  temperature  history  profile,  a  non¬ 
linear,  transient  heat-flow  finite  element  analysis  is  performed  on  the  plane  perpendicular  to  the 
welding  direction. 


Thermal  Analysis 

The  numerical  implementation  of  the  history  dependent  (transient)  heat  transfer  problem  involves 
an  incremental  scheme  with  several  small  time  increments.  This  analysis  procedure  is  addressed  in 
detail  in  references  [27,  28,  29] . 

The  governing  energy  balance  equation  for  transient  heat  transfer  analysis  is  given  as  follows, 
dT 

pCp—(r,  t )  =  -Vr  •  q(r,  t)  +  Q(r,  t)  in  the  entire  volume  Vr  of  the  material  (2.1) 


where  p  is  the  density  of  the  body  ([7,820%/m3]),  Cp  is  the  specific  heat  capacity,  T  is  the 
temperature  ([°C]),  q  is  the  heat  flux  vector,  Q  is  the  internal  heat  generation  rate,  t  is  the  time, 
r  is  the  coordinate  in  the  reference  configuration  and  Vr  is  the  spatial  gradient  operator.  Material 
properties  for  medium  carbon  steel  (A36)  are  used  in  this  study  [30]. 

Convection  boundary  conditions  are  assigned  for  all  free  surfaces.  The  internal  heat  generation 
rate  by  the  welding  torch,  modeled  with  a  ’’double  ellipsoid”  heat  source  model  [16],  is  given  as, 


aba x^/tt 


[f/mm3]  (2,2) 


where  Qb  is  the  welding  heat  input;  rj  is  the  welding  efficiency,  x,  y.  and  2  are  the  local  coordinates  of 
the  double  ellipsoid  model  aligned  with  the  weld  fillet;  a  is  the  weld  width;  b  is  the  weld  penetration; 
c  =  4a  is  the  weld  ellipsoid  length,  and  /  =  0.6  when  the  torch  is  behind  the  analysis  plane,  and 
/  =  1.4  after  the  torch  passes  the  analysis  plane;  v  is  the  torch  travel  speed;  and  t  is  time. 


Elasto-plastic  Mechanical  Analysis 

The  subsequent  history  dependent  stress  analysis  is  performed  by  modelling  the  stress  problem  as  a 
quasi-static  process  in  a  Lagrangian  frame.  This  problem  has  been  covered  by  several  investigators 
[3,  7,  4,  2,  31]. 

The  temperature  values  solved  in  the  previous  thermal  analysis  are  imported  to  the  mechanical 
analysis  as  loading.  Generalized  plane-strain  conditions  are  assumed  to  account  for  the  out-of-plane 
expansion  in  the  structure.  The  longitudinal  (out-of-plane)  strain  ez  is  assumed  to  vary  linearly 
with  x-  and  y-  coordinates  in  the  analysis  plane: 

ez  =  e-x<f>y  +  y<f>x  (2.3) 

where  e  is  the  z-component  of  the  strain  at  the  coordinate  origin  and  the  constants  <j>x  and  <j>y 
represent  the  strain  variations  in  the  y  and  x  axes,  respectively. 
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The  stress  equilibrium  equation  is  given  by, 

Vrcr(r,  t)  +  b(r,  t)  =  0  in  Vr  (2.4) 

where  <r  is  the  stress,  b  the  body  force,  and  t  is  time.  The  mechanical  constitutive  law  is  : 

o'  =  C  (e  —  ep  —  et)  (2.5) 

€P  =  4.a(ff,eq,r)  (2.6) 

f  =  cre  —  <rv  <0  (2.7) 

where  T  is  temperature,  C  is  the  material  stiffness  tensor,  a  is  the  plastic  flow  vector,  e,  ep  and  et 
are  the  total,  plastic  and  thermal  strains  and  eq  is  the  equivalent  plastic  strain.  In  Equation  (2.7), 
/  is  the  yield  function,  ae  is  the  Von  Misses  stress,  and  ay  is  the  yield  stress.  Active  yielding  occurs 
when  f  =  0.  The  mechanical  material  properties  are  assumed  for  A36  [30] . 

Weld  Scaling  Factor 

The  longitudinal  residual  stress  is  positive  at  the  weld  region  and  negative  elsewhere.  This  stress 
distribution  is  caused  by  a  negative  longitudinal  plastic  strain  at  the  weld  region.  Instead  of 
applying  the  exact  plastic  strain  distribution  on  the  3-D  structure,  a  negative  unit  thermal  load  is 
applied  at  the  welding  region. 

A  2-D  linear  analysis  with  a  unit  thermal  load  is  performed  to  determine  the  scaling  factor  to 
the  unit  load  that  produces  equivalent  longitudinal  stress  to  the  welding  residual  stress.  With  the 
acquired  longitudinal  residual  stresses  at  the  free  edges,  are3  in  the  weld  simulation  and  a\  in  the 
unit  thermal  load  analysis  respectively,  the  scaling  factor  is  evaluated  as 

7  =  OVes/eri  (2.8) 


2.2.2  Structural  Analysis 

3-D  Small  Deformation  Analysis; 

The  small  deformation  (linear)  analysis  is  defined  as  follows  [32,  33]: 

Ku  =  f  (2.9) 

where  K  represents  the  linear  stiffness  matrix,  f  represents  the  nodal  load  vector,  and  u  is  the 
nodal  displacement  vector. 

The  linear  deformation  analysis  is  performed  with  applying  a  negative  thermal  load  (T  =  -7). 
Due  to  the  assumption  of  the  linear  relationship  between  load  and  displacement,  the  analysis  is 
accomplished  by  applying  the  (thermal)  load  in  one  step  to  acquire  the  deformation  magnitude. 
3-D  Large  Deformation  Analysis: 

The  large  (non-linear)  deformation  analysis  is  defined  as  solving  the  following  equation  for  u 
[32,  33] 

(  K  +  KG  )u  =  f  (2.10) 

where  K  and  Kg  represent  the  linear  and  non-linear  (stress  stiffening)  stiffness  matrices, 
respectively,  f  represents  the  nodal  load  vector,  and  u  is  the  nodal  displacement  vector. 

A  negative  thermal  load  (T  =  —7)  is  applied.  Unlike  the  linear  analysis  which  can  be 
accomplished  in  one  step,  the  large  deformation  analysis  is  performed  incrementally  with  the  applied 
load  divided  into  small  increments. 
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3-D  Eigenvalue  Buckling  Analysis: 

The  elastic  instability  problem  is  defined  as  an  eigenvalue  problem  as  follows  [32,  33] 

det  (  K  +  AKg  )  =  0  (2.11) 

where  K  and  Kg  are  the  linear  and  non-linear  (stress  stiffening)  stress  stiffness  matrices,  and  X  is 
the  eigenvalue. 

A  3-D  eigenvalue  analysis  is  performed  on  the  structural  model  with  a  unit  negative  thermal 
load  applied  in  the  weld  region  (T  =  -1.0)  to  model  the  uniform  compressive  longitudinal  plastic 
strain  field  occurring  in  welding.  The  eigenvalues  (A*)  represent  the  multipliers  (scaling  factors) 
which  result  in  the  critical  buckling  stress  field  (£rtT)i  when  multiplied  with  the  stress  field  resulting 
from  the  unit  thermal  load 

(crcrX  =  Xi  •  a L  [MPa]  (2.12) 

The  buckling  distortion  is  determined  from  the  eigenvectors  (mode  shapes)  of  the  structure. 
The  structure  may  buckle  in  any  of  the  modes  with  critical  stresses  lower  than  the  residual  stress 
field  due  to  welding.  It  will  tend  to  buckle  with  the  permissible  buckling  mode  having  the  lowest 
critical  stress.  The  permissibility  of  the  modes  is  determined  by  the  constraints  on  the  structure.  If 
certain  buckling  modes  are  suppressed  by  the  mechanical  fixturing  applied,  the  structure  will  tend 
to  buckle  the  next  available  (higher)  mode. 

2.3  Maglev  guideway  beam  modelling 

Dimensions  of  the  Model 

The  actual  length  of  the  guideway  beam  utilized  in  this  project  is  61.92  m.  As  the  beam  has 
a  uniform  cross  section  and  comists  of  alternating  diaphragm  and  crossbeam  stiffeners  at  equally 
spaced  intervals  of  approximately  3  m,  only  a  portion  of  the  beam  is  analyzed  to  simplify  the 
analysis.  The  portion  analyzed  has  a  length  of  9.804m  and  contains  one  end-bearing  diaphragm, 
2  crossbeam  stiffeners  and  one  large  diaphragm.  The  boundary  conditions  used  are  as  follows:  1) 
The  x-  translation  displacement  is  set  to  be  0  on  node  a  (Figure  2.5)  2)  The  x-  and  y-  translation 
displacements  are  set  to  be  0  at  node  b  (Figure  2.5)  3)  The  z-translation  displacement  and  the  x- 
and  y-  rotational  displacements  are  set  to  0  on  all  nodes  at  the  cut  plane  (symmetry  B.C.) 

The  entire  length  of  maglev  guideway  beam,  utilized  in  this  analyzed  project  is  200  feet  (61 
meters)  long.  In  order  to  simplify  the  analysis  model,  a  portion  of  length  9.804  m  which  contains 
one  end-bearing  diaphragm,  one  diaphragm  and  two  cross-beam  stiffeners  inside  the  portion  is 
considered.  In  the  analyzed  models,  the  boundary  conditions  are  as  follows  :  1)  x-translational 
displacement  is  set  to  be  0  at  node  a, (Figure  2.5)  2)  x-  and  y-  translational  displacements  are  set  to 
be  0  at  node  b, (Figure  2.5)  3)  z-translational  displacement  and  x-  and  y-  rotational  displacements 
are  set  0  at  all  nodes  on  the  cut  (symmetry  b.c.)  plane. 

Welding  Sequence 

As  the  current  manufacturing  approach  is  to  tack  weld  the  stiffeners  and  the  web  plates  onto 
the  top  flange  plate  before  welding,  the  same  is  assumed  in  the  analysis.  Since  the  structure  is 
elastic,  the  model  is  analyzed  with  the  static  stresses  applied  simultaneously  to  the  structure. 

Analysis  Cases 

To  examine  the  effect  of  the  guide  rails  on  the  maglev  guideway  beam  if  the  structure  is 
manufactured  sequentially,  analyses  are  performed  for  two  models:  one  with  and  the  other  without 
guide  rails.  The  analysis  models  are  shown  in  Figures  2.5  and  2,6,  respectively. 
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2.3.1  2-D  model  and  welding  simulation  Simulation 

The  2-D  thermo-mechanical  analysis  is  performed  to  determine  the  scaling  factor  from  welding  to 
be  applied.  The  boundary  conditions  utilized  in  this  analysis  are  shown  in  Figure  2.3.  The  inclined 
plate  represents  the  12  mm  web  plate  and  the  bottom  plate  represents  the  30mm  bottom  flange. 
The  heat  from  welding  at  the  fillet  is  applied  as  a  heat  source  as  per  Equation  (2.2),  All  the  free 
surfaces  are  the  convective  surfaces.  The  node  on  the  bottom  left  comer  is  restrained  in  all  degrees 
of  freedom  and  the  node  on  the  web  plate  top  canter  is  restrained  in  the  x  direction  to  prevent 
rotation.  The  2-D  finite  element  mesh  in  the  heat  transfer  and  mechanical  analysis  of  the  panel  is 
illustrated  in  Figure  2.4.  The  model  consists  of  339  8-noded  quadrilateral  quadratic  elements  with 
1136  nodes.  The  finite  element  solutions  are  performed  using  an  SMP  in-house  Fortran90  finite 
element  code. 

Weld  Simulation 


Figure  2.3:  Boundary  Conditions  for  2-D  welding  simulation  models 
Thermo-mechanical  analysis  of  unit  thermal  load 

2-D  generalized  plane  strain  models  are  developed  to  compute  the  stress  <tl  resulting  from  a 
unit  thermal  load.  The  same  2D  mesh  as  in  earlier  case  is  used  for  this  analysis.  A  negative  unit, 
load  is  applied  in  the  weld  zone.  Ambient  temperature  is  applied  at  all  other  nodes. 

2.3.2  3-D  Structural  Analysis 

Two  cases  were  analyzed:  The  first  model  had  the  guide  rails  and  consisted  of  a  total  of  1691  nodes 
and  1948  elements,  including  1658  (4-noded)  shell  elements  and  290  truss  elements  (Figure  2.5).  The 
second  model  did  not  have  the  guide  rails  and  consisted  of  a  total  of  1539  nodes  and  1726  elements, 
including  1510  shell  elements  and  216  truss  elements  (Figure  2.6).  The  truss  elements  represent 
the  weld  zone  and  the  thermal  loads  are  applied  to  the  truss  elements.  Both  the  deformation  and 
buckling  are  performed  on  both  models  using  the  ABAQUS  finite  element  software. 
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Figure  2.4:  FEA  mesh  for  2-D  welding  process  analysis 

The  thermal  loads  are  applied  at  the  truss  elements.  Based  on  the  weld  size,  the  cross-section 
area  of  the  truss  elements  are  32mm2  along  the  line  connecting  guide  rails  and  the  top  plate,  18mm2 
along  other  welding  lines  except  the  elements  at  the  cut  (symmetry  B.C.)  plane  which  has  9mm2, 
The  deformation  and  buckling  analyses  are  performed  by  using  ABAQUS  software.  No  preload  is 
applied  in  the  analyses. 


2.4  Results 

2.4.1  2-D  Analyses 

Weld  Simulation 

Figure  2,7  shows  the  distribution  of  residual  stresses  in  Z  direction  as  computed  by  the  2-D 
simulation  of  welding.  The  stresses  are  tensile  in  the  weld  zone  and  in  the  heat  affected  zone. 
However,  away  from  the  weld  zone,  the  stress  changes  the  sign  from  tensile  to  compressive.  The 
peak  residual  longitudinal  stress  at  the  weld  center  is  412  MPa  and  residual  stress  at  the  free  edge 
is  (Tres  =  -12  MPa.  Figure  2.8  shows  the  plastic  strain.  All  of  the  plastic  strains  are  negative. 
Unit  Load 

Figure  2.9  shows  the  distribution  of  stresses  in  Z  direction  for  the  2D  thermo-mechanical  analysis 
due  to  the  unit  load.  The  maximum  residual  stress  at  the  weld  center  is  0.94  MPa  and  that  at  the 
free  edges  is  a\  =  -0.0322  MPa.  Thus,  the  multiplier  7  (scaling  factor)  for  welding  applied  on  the 
model  can  be  acquired  as 


7  = 


&res 

G\ 


-12 

-0.0322 


=  372 


(2.13) 
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2.4.2  3-D  Analysis 

Model  with  Guide  Rails 
Deformation  Analysis 

Figure  2.10  illustrates  the  small  deformation  and  large  deformation  analysis  results  with  20X 
magnification.  The  linear  analysis  captures  the  bowing  effect  at  the  free  end  side.  The  large 
deformation  analysis  captures  the  combined  effects  of  bowing  and  buckling. 

Significant  distortion  is  observed  at  the  free  end.  Figure  2.11  shows  the  x-direction  displacements 
along  the'  guide  rail  edges,  i.e.  along  line  A  and  line  B  shown  in  Figure  2.10.  The  y-direction 
displacements  along  both  longitudinal  lines  are  shown  in  Figure  2.12. 

Buckling  Modes 

The  first  four  buckling  modes  of  the  model  with  guide  rails  are  shown  in  Figure  2.13,  where  it  is 
observed  that  the  buckling  deformation  on  the  web  flange  are  more  evident  than  on  the  top  flange. 

The  first  and  second  eigenmodes  can  be  categorized  as  one  group  with  repeated  eigenvalues, 
because  those  are  basically  the  same  type  of  buckling  waves.  The  major  difference,  as  shown  in 
Figure  2.14,  is  that  the  first  mode  is  symmetric  about  the  YZ-plane  while  the  second  mode  is 
anti-symmetric. 
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Figure  2,6:  FEA  mesh  for  3-D  structural  analysis  in  the  case  of  without  guide  rails 


Model  without  Guide  Rails 
Deformation  Analysis 

The  deformations  in  the  model  without  guide  rails  are  shown  in  Figure  2.15  for  large  deformation 
and  linear  analysis,  respectively.  Unlike  the  results  in  the  case  with  guide  rails,  the  deformation 
result  of  linear  analysis  shows  that  there  is  no  bowing  effect,  because  no  welding  load  is  applied  on 
the  lines  designed  to  be  connected  with  guide  rails. 

Buckling  waves  near  the  edges  are  observed  on  the  top  panel  in  the  large  deformation  analysis 
result,  because  there  are  no  constraints  applied  by  the  guide  rails.  Line  C  in  Figure  2.15  represents 
the  line  where  the  guide  rails  are  to  be  welded.  Figure  2,16  shows  the  x-  and  y-  direction 
displacements  along  the  edge  of  the  top  plate,  along  the  Line  C  in  Figure  2,15. 

Buckling  Modes 

The  first  four  buckling  modes  of  the  model  with  guide  rails  are  shown  in  Figure  2.17.  Similar 
to  the  modes  in  the  analysis  case  with  guide  rails,  the  buckling  deformations  are  more  evident  on 
the  side  panels  than  on  the  top  panel.  However,  compared  with  Figure  2.13,  because  guide  rails 
are  not  present,  the  top  panel  has  more  deformation. 

Similarl  to  the  model  with  guide  rails,  repeated  eigenvalues  and  corresponding  eigenmodes  are 
also  observed. 
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Mode 

w/  Guide  Rails 

w/o  Guide  Rails 

Mode 

w/  Guide  Rails 

w/o  Guide  Rails 

i 

98.989" 

114.92 

16 

159.34 

181.92 

100.09 

116.53 

161.44 

182.96 

108.00 

124.01 

165.70 

189.55 

108.79 

125.03 

191.35 

217.68 

110.99 

127.04 

194,15 

221.67 

6 

111.99 

128.05 

21 

222.18 

255.87 

113.92 

129.78 

223.15 

257.44 

116.75 

132.74 

234.63 

269.67 

116.97 

133.40 

235.69 

271.14 

119.77 

138.63 

266.78 

278.04 

11 

124.45 

143.79 

26 

269.87 

297.88 

125,61 

145.91 

270.71 

298.87 

150.75 

168.38 

274.10 

307.73 

151.89 

173.36 

282.12 

320.78 

156.32 

174.18 

295,09 

327.11 

Table  2.1:  Comparison  of  eigenvalues  in  buckling  eigenvalues  below  the  scaling  factor  between 
different  models 


Mode 

w /  Guide  Rails 

w/o  Guide  Rails 

31 

295.95 

331.89 

298.86 

337.94 

305.37 

342.94 

306.03 

351.36 

307,56 

354.19 

36 

313.16 

362.72 

317.88 

368.26 

318.29 

370.79 

320.21 

371.83 

321.75 

- 

41 

332.54 

- 

334.63 

- 

340.48 

347.23 

348.47 

- 

46 

350.51 

351.56 

- 

355.92 

- 

360.37 

- 

361.87 

- 

51 

367.80 

- 

Table  1:  Comparison  of  eigenvalues  in  buckling  eigenvalues  below  the  scaling  factor  (Cont’d) 
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Figure  2.8:  Plastic  Strain  obtained  from  2-D  analysis 
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Figure  2.9:  Longitudinal  Residual  Stress  using  Unit  Load  obtained  from  2-D  analysis  (MPa) 
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Figure  2.12:  Y-direction  Displacements  of  the  Guide  Rails 
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Symmetry  B,C. 


Symmetry  B,C, 


The  first  buckling  mode  with  top  view 


The  second  buckling  mode  with  top  view 


Figure  2,14:  First  and  Second  Buckling  Modes  with  Top  View 
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Symmetry  &C,  Symmetry  B.C. 


Second  buckling  mode  (X  =  1 1 6.53)  Fourth  buckling  mode  (X  - 1 25.03) 


Figure  2.17;  The  First  Four  Buckling  Modes  for  the  Model  without  Guide  Rails 


Chapter  3 


Domain  Decomposition  and 
Adaptivity  Methods 
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3.1  Domain  Decomposition  Method  for  Parallel  Computing  and 
Sensitivity  Analysis 


Motivation  of  Domain  Decomposition  Method: 


•  The  practical  engineering  simulation  problem  is  very  large  and  complex 

•  The  single  computer  power  is  limited  and  can  not  handle  such  a  big  job 

•  Parallel  processing  is  introduced  to  resolve  this  issue 

•  Domain  decomposition  method  is  a  native  way  to  split  problems  for  parallel  processing 


The  nature  of  Dd  method  is  divide  and  conquer.  The  domain  is  decomposed  into  several 
subdomains  and  they  are  assigned  to  individual  processor.  The  shared  boundary,  or  its  related 
Lagrange  multiplier  equation  is  solved  first,  then  the  internal  nodes  in  each  subdomain.  FETI 
methods  is  a  family  of  iterative  domain  decomposition  methods,  it  includes: 


•  FETI  method  [34] 

•  FETI-2  method  [35,  36] 

•  FETI-Dual  Primal  method  [37,  38] 


The  FETI  Family  Methods  are  numerically  scalable  for  second  and  fourth  order  elasticity 
problems,  such  as  plane  stress,  plane  strain,  solid  mechanics,  plate  and  shell  problems.  Each 
subdomain  contains  the  nodes  from  the  shared  boundary,  so  connectivity  between  neighboring 
domains  are  enforced  by  introducing  Lagrange  multiplier,  which  forms  the  basic  unknowns  of 
interface  equation.  In  FETI-Dual  Primal  methods,  the  displacement  and  rotation  degrees  of 
freedom  from  some  corner  nodes  is  also  introduced  as  unknowns.  Finally,  the  above  equations 
is  solved  by  an  iterative  method,  such  as  preconditioned  conjugate  projected  gradient  method 
(PCPG). 

When  doing  sensitivity  analysis,  since  only  the  right  hand  side  term  is  changed,  it  is  preferred 
to  have  the  direct  inverse  of  stiffness  matrix.  FETI  methods  is  not  suitable  here  since  it  is  based 
on  iterative  method.  A  similar  domain  decomposition  method  is  used  to  make  the  calculation  of 
inverse  matrix  easier,  now  all  the  degrees  of  freedom  from  the  shared  boundary  forms  the  basic 
unknowns,  and  the  corresponding  equation  is  called  Schur  equation.  The  inverse  of  global  stiffness 
matrix  can  be  replaced  by  the  inverse  of  each  subdomain  stiffness  matrix  and  the  inverse  of  Schur 
matrix. 

An  example  of  two  subdomains  is  illustrated  in  Figure  3.1 


Subdomain  1 


Subdomain  2 


The  finite  element  equation  is  as  follows: 

■All  0  A  i3  ui  /i 

0  Ag  2  A23  u2  =  /2 

.  A13  A23  -4  3 3  fj,  fB 

The  Schur  equation  is  derived  by  performing  Gauss  elimination  on  ui,u2 


Sli  -  /b 


where, 


3  —  A33  .4 13  —  A23  A221  A23 

f  B  =  jb  —  Af3rin1  fi  —  A23A221  fi 
Once  fj,  is  solved,  we  can  derive  ui  ,«2  from  equation  (3.1) 

ui  =  An  (/i  -  A13M)  ^2  =  A22H/2  -  A23ju) 

To  obtain  the  inverse  matrix  of  this  finite  element  equation,  which  is  the  same  as  solving  this 
equation,  we  need  to  calculate  three  inverse  matrices 

AJi1 ,  A22 ,  and  S_1 
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this  can  be  done  by  most  direct  solvers,  like  Watson  solver  and  ESSL  solver 
Future  Work  on  this  topic  will  include  the  following: 

•  Testing  on  large  and  very  large  models 

•  Programming  for  sensitivity  analysis 

•  Implement  on  distributed  machines  by  using  Message  Passing  Interface  (MPI) 


3.2  Adaptivity  for  Welding  Simulations 

Research  has  started  for  the  implementation  of  adaptivity  in  welding  simulations  of  large  structures. 
A  flow  chard  of  the  approach  is  illustrate  in  Figure  3.2.  The  approach  is  based  on  h-adaptivity 
using  Lagrange  multipliers  to  constrain  incompatible  elements. 

(£t)("Hfo)  (-) 

where,  the  original  equation  when  there  was  no  h-adaptivity  applied  is  K  *  u  =  f ,  A  is  the 
Lagrange  multipliers.  The  operator  A  is  assembled  with  columns  based  on  the  mathematical 
relations  of  the  dependent  nodes  among  other  nodes.  For  example,  if  unodel  =  |(unode2  +  unode3), 

1  ...  \ 

-0.5  ... 

A==  -0.5  ... 

\  i  -.) 

After  rearrangement,  Equation  (3.3)  will  becomes 

(  0  -A^A  )  (  A  )  =  (  -AtK ~H>  )  <3*4> 

the  left-hand  side  is  in  a  form  of  an  upper-triangle  matrix.  With  the  A  being  evaluated,  the  u 
is  then  solved  by  / 


u  =  K_1(f/  —  AA)  (3.5) 

Thus,  the  equation  for  applying  h-adaptivity  can  be  solved  in  an  efficient  way. 

Further  work  is  to  test  the  whole  procedure,  combining  the  Lagrange  multiplier  procedure  and 
refining  elements,  with  heat  transfer  models  and  mechanical  models. 
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1  Technical  Objectives 


The  objective  of  this  project  is  to  develop  a  computer  aided  weld  design  tool  for  a  computer 
integrated  design  and  manufacturing  system. 

The  weld  design  tool  will  allow  for  the  determination  of  appropriate  welding  conditions  and 
if  needed  auxiliary  heating  for  a  specific  structural  design  such  that  welding  distortion  is  within 
tolerance.  The  database  will  be  integrated  with  Computer  Aided  Production  Engineering  (CAPE) 
software. 

The  project  contributes  to  a  collaborative  effort  by  Maglev  Inc  and  other  Universities  entitled 
’’Demonstration  of  Computer  Aided  Manufactuting  Techniques  for  the  Precision  Fabrication  of 
Large  Steel  Curved  Plate  Beam  Components  for  Shipbuilding  and  Other  Industries.” 


2  Technical  Approach 


The  proposed  program  is  organized  in  the  following  five  tasks: 

1.  Identification  of  structural  features  that  welding  distortion 

2.  Inverse  method  and  database  architecture 

3.  Generation  of  a  database 

4.  Integration  with  computer  aided  production  engineering  software 

5.  Modifications  and  refinement  using  production  tests 


3  Progress 

3.1  Identification  of  structural  features  that  welding  distortion 


Welding,  among  all  mechanical  joining  processes,  has  been  employed  at  an  increasing  rate  for 
its  advantages  in  design  flexibility,  cost  savings,  reduced  overall  weight  and  enhanced  structural 
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performance.  However,  welding  induces  various  types  of  distortions  as  discussed  in  detail  by 
Masubuchi  [1],  To  access  the  effects  of  welding  on  structure  efficiently,  and  hence  in  turn  to 
implement  various  distortion  mitigation  techniques,  a  validated  method  for  predicting  welding 
induced  distortion  is  necessary. 

Thinner  section  components  made  of  higher  strength  steels  are  being  commonly  utilized  in 
shipbuilding,  railroad  and  aerospace  industries  in  fabricating  large  structures  to  achieve  reduction  in 
overall  weight  and  more  controllable  manufacturing.  However,  for  structures  made  of  relatively  thin 
components,  welding  can  introduce  significant  buckling  distortion  which  causes  loss  of  dimensional 
control,  structural  integrity  and  increased  fabrication  costs  due  to  poor  fit-up  between  panels.  A 
predictive  analysis  technique  can  determine  the  susceptibility  of  a  particular  design  to  buckling 
distortion.  Further,  a  predictive  analysis  tool  can  assist  in  the  selection  of  geometry  and  welding 
conditions  that  will  induce  minimum  distortion.  Flame  straightening  is  the  commonly  used 
technique  to  correct  the  out-of-plane  distortion  resulting  from  welding  processes,  but  this  is  an 
inaccurate,  labor  intensive  and  costly  process.  Also  it  is  rather  a  corrective  action  after  the  damage 
is  done  and  not  a  preventive  measure  which  is  generally  desired  in  engineering  processes. 

Finite  element  techniques  have  been  used  in  the  prediction  of  welding  residual  stress  and 
distortion  for  more  than  two  decades.  Due  to  the  nature  of  the  process,  additional  complexities  are 
involved  in  the  FEA  of  welding  compared  to  traditional  mechanics,  such  as  temperature  and  history 
dependent  material  properties;  high  gradients  of  temperature,  stress  and  strain  fields  with  respect  to 
both  time  and  spatial  coordinates;  large  deformations  in  thin  structures  and  phase  transformation 
and  creep  phenomena. 

Earlier  studies  of  welding  accounted  for  the  non-linearities  due  to  temperature  dependent 
material  properties  and  plastic  deformations  [2,  3,  4],  The  majority  of  those  analyses  were  limited 
to  two-dimensions  on  the  plane  perpendicular  to  the  welding  direction,  but  good  correlations 
have  been  observed  between  the  numerical  predictions  and  experimental  results  [5,  6,  7,  8],  and 
especially  for  residual  stress  predictions,  2-D  models  provided  accurate  estimations  comparable  to 
3-D  analyses,  since  the  stress  field  exhibits  a  fairly  uniform  distribution  through  the  length  of  the 
work-piece.  Argyris  et.  al.  [5]  computed  the  thermo-mechanical  response  using  2-D  models  in  a 
staggered  solution  strategy  to  combine  and  integrate  the  thermal  and  mechanical  computational 
steps.  Rybicki  et.  al.  [6]  performed  thermo-elasto-plastic  analysis  on  a  2-D  axisymmetric  finite 
element  model  for  a  two-pass  girth-butt  welded  pipe  problem,  and  verified  the  numerical  results 
with  the  experimentally  obtained  temperature  history  and  residual  stress  distributions.  Papazoglu 
and  Masubuchi  [7]  solved  the  multipass  GMAW  process  problem  by  performing  uncoupled  2-D  heat 
transfer  and  stress-strain  analyses,  incorporating  the  phase  transformation  strains. 

2-D  models,  as  mentioned  above,  have  been  particularly  useful  with  their  high  efficiency  and 
accuracy  in  determining  the  solution  in  the  analysis  plane  and  reduced  computational  requirements. 
However,  for  welding  practices  where  tack  welding  or  fixturing  allow  out-of-plane  movement  2-D 
analyses  may  not  be  accurate,  particularly,  in  distortion  predictions  [9].  Furthermore,  longitudinal 
heat  transfer  and  instability  aspects,  and  end  effects  (i.e.  due  to  initiation  and  termination  of  the 
heat  source)  cannot  be  realized  in  two  dimensional  formulations. 

Most  of  the  currently  performed  welding  simulations,  both  2-D  and  3-D,  are  based  on  small 
deformation  assumption  and  are  limited  to  simpler  structures  and  weld  geometries  (e.g.  butt  joints) 
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or  focusing  only  to  the  heat  affected  zones,  ignoring  the  surrounding  structure.  A  small  deformation 
analysis  assumes  infinitesimal  displacements  and  loads  being  applied  to  the  undeformed  geometry. 
The  interaction  between  the  weld  zone  and  the  structure  is  effective  on  the  accumulated  distortion, 
and  large  deformation  modes  in  unrestrained  structures  may  not  be  captured  with  this  type  of 
analysis  [9],  [10].  Brown  and  Song  [9]  have  performed  2-D  axisymmetric  and  3-D  weld  simulations 
of  a  ring  stiffened  cylinder  structure,  and  concluded  that  2-D  analysis  overestimated  the  rotation  of 
the  ring  dining  the  heating  segment,  and  it  was  very  sensitive  to  model  modifications,  such  as  joint 
clearance  and  location  of  constraints.  Michaleris  et.  al.  [10]  studied  the  effects  of  the  restraints 
and  the  solidified  portions  of  the  weld  on  the  residual  stress  and  distortion  profiles  by  comparing 
the  performance  of  2-D  and  3-D  weld  simulations. 

Oddy  et.  al.  [11]  examined  the  butt  welding  of  a  bar  via  3-D  FEM,  and  computed  the 
temperature,  strain  and  stress  fields,  Tekriwal  and  Mazumder  [12,  13]  simulated  thermal  and 
elasto-plastic  response  of  the  butt-welded  plates  through  3-D  models,  considering  filler  material 
addition.  Multi-pass  welding  simulation  of  plates  and  experimental  validation  have  been  addressed 
in  [14,  15]. 

Welding-induced  buckling  of  thin-walled  structures  has  been  investigated  in  greater  detail  by 
[16,  17,  18].  Ueda  et.  al.  [16,  19]  presented  a  methodology  to  determine  the  buckling  behavior  of 
plates  by  large  deformation  elastic  FEA  and  employing  inherent  strain  distributions.  Michaleris  et. 
al.  [17,  20]  developed  a  predictive  buckling  analysis  technique  for  thin  section  panels,  combining 
decoupled  weld  process  simulations  and  eigenvalue  buckling  analyses.  Tsai  et.  al.[18j  studied  the 
distortion  mechanisms  and  the  effect  of  welding  sequence  on  panel  distortion. 

For  the  welding  practices  where  tack  welds  or  fixturing  are  used  to  restrict  the  movement  of  the 
welded  parts,  the  structural  response  may  be  evaluated  by  means  of  decoupled  2-D  welding  and 
3-D  buckling  simulations.  When  mechanical  fixturing  on  the  structure  prevents  the  longitudinal 
shrinkage  during  welding,  the  out-of-plane  structural  behavior  doesn’t  have  influence  on  the  in¬ 
plane  welding  response,  and  buckling  is  only  observed  after  the  restraints  are  removed  and  the 
structure  cools  down.  Exploiting  this  fact,  Michaleris  et.  al.  [17]  proposed  the  aforementioned 
buckling  prediction  technique  with  uncoupled  weld  simulation  and  structural  buckling  analyses. 
They  expressed  the  residual  stress  profile  from  the  2-D  welding  simulations  as  buckling  stress  on 
the  3-D  structural  model.  This  approach  is  analogous  to  the  work  by  Ueda  et.  al.  [19],  where  the 
concept  of  inherent  strain  is  used  to  generate  the  welding  residual  stresses  by  applying  a  prescribed 
thermal  strain  field  using  empirical  methods.  In  the  former  study,  however,  residual  stresses  are 
calculated  with  weld  process  simulations,  which  provides  improved  estimations  for  buckling  analysis 
compared  to  empirically  determining  the  residual  stress. 

Phase  transformations  and  transformation  plasticity  have  also  been  incorporated  in  the  analysis 
as  recent  developments  [11,  21,  22].  The  primary  objective  there  is  to  more  accurately  model  the 
residual  stress  distribution,  microstructure  and  local  distortion  in  the  area  immediately  adjacent 
to  the  weld. 

In  this  work  the  decoupled  2-D  and  3-D  finite  element  analysis  technique  by  Michaleris  et. 
al.  [20,  17]  is  applied  to  evaluate  welding-induced  buckling  of  the  welded  panels.  Effects  of  the 
following  process  and  design  parameters  are  investigated: 
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•  panel  cross-section  (panel  thickness) 

•  panel  size, 

•  welding  heat  input, 


3.1.1  Analysis  Approach  -  Modelling  the  Welding  Distortion 


Following  the  work  of  Michaleris  et.  al.  [20,  17],  the  response  of  welded  panels  is  evaluated  in  two 
steps  by  combining  two-dimensional  welding  simulations  with  three-dimensional  structural  analyses 
in  a  decoupled  approach. 


2-D  Thermo-mechanical  Weld  Simulation  : 


A  two  dimensional  thermo-elasto-plastic  analysis  is  performed  to  determine  the  angular 
distortions,  residual  stresses,  and  plastic  strain  fields  during  the  welding  process  ignoring  the 
structural  response.  Residual  stresses  are  caused  by  the  negative  plastic  strains  resulting  from 
the  welding  thermal  cycle. 


3-D  Eigenvalue  Analysis  : 


The  buckling  distortion  and  critical  buckling  stresses  are  consequently  determined  by  an 
eigenvalue  analysis  applying  the,  mostly  uniform  and  compressive,  longitudinal  plastic  strain  field 
of  the  2-D  weld  model  on  the  3-D  structural  model  as  equivalent  load. 

A  constant,  negative  thermal  load  is  applied  at  the  weld  region  to  introduce  the  effects  of 
welding  into  the  3-D  structure.  Thermal  loading  is  used  rather  than  mapping  the  plastic  strain 
field,  which  would  require  a  complex  analysis  procedure.  An  eigenvalue  analysis  is  performed  to 
determine  the  critical  residual  stresses  and  buckling  distortions. 


Welding  Simulation 


The  welding  simulation  involves  a  thermal  and  a  mechanical  analysis.  The  effect  of  mprVianjr-ai 
response  is  assumed  to  be  negligible  on  the  thermal  behavior,  thus  the  temperature  field  is  solved 
independently  from  the  mechanical  solution.  To  determine  the  temperature  history  profile,  a  non¬ 
linear,  transient  heat-flow  finite  element  analysis  is  performed  on  the  plane  perpendicular  to  the 
welding  direction. 


Thermal  Analysis 
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The  numerical  implementation  of  the  history  dependent  (transient)  heat  transfer  problem 
involves  an  incremental  scheme  with  several  small  time  increments.  The  solution  at  a  given  time 
increment  is  obtained  by  using  the  solution  at  the  previous  time  increment  as  an  initial  condition. 
This  problem  is  addressed  in  detail  in  references  [3,  23,  12]. 

The  governing  energy  balance  equation  for  transient  heat  transfer  analysis  is  given  as  follows, 
dT 

pCp^-(r,  t)  =  -Vr  •  q(r,  t )  +  Q(r,  t)  in  the  entire  volume  Vr  of  the  material  (1) 

where  p  is  the  density  of  the  body  ([7,82 Okg/m3]),  Cp  is  the  specific  heat  capacity,  T  is  the 
temperature  ([°C]),  q  is  the  heat  flux  vector,  Q  is  the  internal  heat  generation  rate,  t  is  the  time, 
r  is  the  coordinate  in  the  reference  configuration  and  Vr  is  the  spatial  gradient  operator.  Material 
properties  for  medium  carbon  steel  (A36)  are  used  in  this  study  (Figure  1). 


Figure  1:  Conductivity  (k),  specific  heat  (Cp),  and  air  convection  (h)  for  A36. 

The  nonlinear  isotropic  Fourier  heat  flux  constitutive  relation  is  enforced;  using  the  temperature- 
dependent  thermal  conductivity,  k  ([og^^a]), 

q  =  -kVrT  {W/mm2}  (2) 


Convection  boundary  conditions  are  assigned  for  all  free  surfaces.  The  internal  heat  generation 
rate  by  the  welding  torch,  modeled  with  a  ’’double  ellipsoid”  heat  source  model  [24],  is  given  as, 


abaK^K 


[W/mm3]  (3) 
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where  Qb  is  the  welding  heat  input;  rj  is  the  welding  efficiency,  x,  y.  and  2  are  the  local  coordinates  of 
the  double  ellipsoid  model  aligned  with  the  weld  fillet;  a  is  the  weld  width;  b  is  the  weld  penetration; 
c  =  4a  is  the  weld  ellipsoid  length,  and  /  =  0.6  when  the  torch  is  behind  the  analysis  plane,  and 
/  =  1-4  after  the  torch  passes  the  analysis  plane;  v  is  the  torch  travel  speed;  and  t  is  time. 


Elasto-plastic  Mechanical  Analysis 


The  subsequent  history  dependent  stress  analysis  is  performed  by  modelling  the  stress  problem 
as  a  quasi-static  process  in  a  Lagrangian  frame.  This  problem  has  been  covered  by  several 
investigators  [5,  15,  4,  13,  11).  Similar  to  the  heat  transfer  analysis,  the  numerical  implementation 
of  the  quasi-static  analysis  involves  an  incremental  scheme  with  several  small  static  increments.  The 
solution  at  a  given  time  interval  is  obtained  by  using  the  solution  at  the  previous  time  increment 
as  an  initial  condition. 

The  temperature  values  solved  for  in  the  previous  thermal  analysis  are  imported  to  the 
mechanical  analysis  as  loading.  Generalized  plane-strain  conditions  are  assumed  to  account  for 
the  out-of-plane  expansion  in  the  structure.  The  longitudinal  (out-of-plane)  strain  is  assumed  to 
vary  linearly  with  x-  and  y-  coordinates  in  the  analysis  plane: 

€2  e  X(f)y  -f-  yd>x  (4) 

where  e  is  the  z-component  of  the  strain  at  the  coordinate  origin  and  the  constants  <f>x  and  <py 
represent  the  strain  variations  in  the  y  and  x  axes,  respectively. 

The  stress  equilibrium  equation  is  given  by, 

Vro(r,  t)  +  b(r,  t)  =  0  in  Vr  (5) 

where  o  is  the  stress,  b  the  body  force,  and  t  is  time.  The  mechanical  constitutive  law  is  : 

&  =  c  (e-ep-it)  (6) 

dp  =  dq  •  a  (tr,  eq,  T)  (7) 

/  =  ae  -  oy  <  0  (8) 

where  T  is  temperature,  C  is  the  material  stiffness  tensor,  a  is  the  plastic  flow  vector,  e,  ep  and  et 
are  the  total,  plastic  and  thermal  strains  and  eq  is  the  equivalent  plastic  strain.  In  Equation  (8),  / 
is  the  yield  function,  ae  is  the  Von  Misses  stress,  and  ay  is  the  yield  stress.  Active  yielding  occurs 
when  /  =  0.  Figures  2  and  3  illustrated  the  mechanical  material  properties  assumed  for  A36. 


Structural  Analysis 


The  longitudinal  residual  stress  distribution  (oy)  computed  in  the  2-D  analyses  are  compared 
to  the  critical  buckling  stresses  (er CT)  of  the  structure  from  the  3-D  structural  analysis  to  determine 
if  the  structure  will  buckle. 
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Figure  5:  FEA  mesh  for  3-D  buckling  analysis 
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The  structural  analysis  involves  elastic  eigenvalue  analysis.  Incremental  large  deformation 
analyses  may  also  be  performed  to  determine  the  onset  of  buckling,  buckling  and  post  buckling 
stages  in  response  to  increasing  stress,  but  they  are  computationally  intensive  and  are  usually  used 
for  validating  the  predictive  methodology  [17]. 


3-D  Eigenvalue  Buckling  Analysis: 


The  elastic  instability  problem  is  defined  as  an  eigenvalue  problem  as  follows 

det  (  K  +  AKq  )  =  0  (9) 

where  K  and  KG  are  the  linear  and  non-linear  strain  stiffness  matrices,  and  A  is  the  eigenvalue. 

A  3-D  eigenvalue  analysis  is  performed  on  the  structural  model  with  a  unit  negative  thermal 
load  applied  in  the  weld  region  (T  =  -1.0)  to  model  the  uniform  compressive  longitudinal  plastic 
strain  field  occurring  in  welding.  The  eigenvalues  (A *)  represent  the  multipliers  (scaling  factors) 
which  result  in  the  critical  buckling  stress  field  (o-cr)i  when  multiplied  with  the  stress  field  resulting 
from  the  unit  thermal  load  (ctl).  Equation  (10)  shows  the  computation  of  the  critical  residual  stress 
distribution  at  the  plate  midspan. 

(o’cr)i  =  Aj  •  cri  [MPa]  (10) 


The  buckling  analyses  may  yield  negative  eigenvalues,  which  often  cannot  be  explained  by 
!!  physical”  behavior.  Those  situations  can  be  avoided  by  applying  enough  preload  Tp,  to  load  the 
structure  just  below  the  buckling  load  before  performing  the  eigenvalue  analysis.  In  such  a  case, 
the  critical  buckling  stress  in  Equation  (10)  is  determined  as 

(o'er)*  =  (  Af  +  |  Tp  |  )  •  a i  [MPa]  (11) 


Buckling  distortion  is  determined  from  the  eigenvectors  (mode  shapes)  of  the  structure.  The 
structure  may  buckle  in  any  of  the  modes  with  critical  stresses  lower  than  the  residual  stress  field 
due  to  welding.  It  will  prefer  to  buckle  with  the  permissible  buckling  mode  having  the  lowest 
critical  stress.  The  permissibility  of  the  modes  are  determined  by  the  constraints  on  the  structure. 
If  certain  buckling  modes  are  suppressed  by  the  mechanical  fixturing  applied,  the  structure  will  tend 
to  buckle  the  next  available  (higher)  mode.  The  weight  of  the  structure  might  have  an  influence  of 
causing  even  higher  buckling  modes. 


3.1,2  Verification 


Experimental  Setup 
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The  predictive  buckling  analysis  approach  taken  in  this  work  is  verified  with  the  experiments. 
The  panels  of  the  same  size  as  modelled  are  used.  All  other  welding  conditions  like  heat  input  and 
boundary  conditions  are  also  used  same  as  in  the  Finite  Element  Model. 

Three  different  panels  with  appropriate  welding  conditions  are  used  in  this  experiment  as  listed 
in  Table  1.  The  welding  setup  consists  of  two  welding  guns  on  either  side  of  the  stiffener  plate.  The 
guns  are  3.5”  offset  from  each  other  with  one  gun  following  the  other.  This  is  shown  in  the  Figure 
6.  Constant  Voltage  Metal  Inert  Gas  welding  is  carried  out.  Welding  electrode  of  0.045”  Dia.  is 
used  with  a  mixture  of  75/25  Argon  and  Carbon  Dioxide  shielding  gas.  The  welding  conditions  are 
set  to  give  short  circuiting  transfer  mode  of  the  welding. 


Figure  6:  Experimental  Arrangement 

A  linear  motion  device  is  used  for  getting  constant  velocity  of  welding  and  also  to  maintain  a 
constant  stickout  through  the  weld.  The  welding  gun  speed  for  all  welds  is  maintained  at  15  Inches 
per  minute  and  an  electrical  stickout  of  1/2”.  The  gas  flow  of  35  cubic  Feet  per  hour  is  used  for 
shielding.  Drag  angle  of  15  degrees  is  used  between  the  gun  and  the  vertical  plane. 

The  stiffeners  are  tack  welded  to  the  base  plate  to  make  the  T-shaped  panels  ready  for  the 
welding.  The  vertical  stiffener  plate  is  held  horizontally  with  screws  and  the  rammed  down  with 
pneumatic  cylinders.  The  panel  is  supported  by  a  copper  plate  during  welding  for  rapid  post¬ 
welding  cooling  and  also  to  avoid  welding  of  specimen  to  the  base  iron  table.  The  edges  of  the 
specimen  are  let  loose  and  they  are  not  restricted  in  any  degree  of  freedom.  Welding  is  carried  out 
on  all  the  panels  shown  in  Figures  7  and  8  with  the  respective  welding  conditions  listed  in  Table  1. 
No  water  or  forced  air  cooling  of  panel  is  done  after  the  welding  is  complete.  The  panel  is  taken 
out  once  it  Is  cooled  down  and  the  distortion  is  examined. 

The  panel  is  fabricated  by  joining  a  small  stiffener  plate  to  a  large  plate  longitudinally  in  a 
T-joint  configuration.  Three  different  panel  geometries  are  considered,  as  illustrated  in  Table  1. 
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Figure  7:  10”  X  10”  Panel  in  1/16”  and  1/8”  thickness 


Case 

Size 

in. 

Thickness 

in. 

Volts 

Amps 

Travel  Speed 
in. /min 

Wire  feed 
in./min 

1 

10x10 

17.8 

15 

75 

2 

10x10 

1/8 

Ell 

167 

15 

180 

3 

12x12 

00 

rHi 

167 

15 

Table  1:  Weld  Parameters 


The  fillet  welds  of  size  1/5  in.  and  1/16  in.  on  thick  and  thin  plates  are  performed  respectively 
on  both  sides  with  dual-torches  by  gas  metal  arc  (GMAW).  The  effect  of  geometry  and  process 
parameters  are  studied. 


3.1.3  Numerical  Implementation 


Figure  9:  Boundary  Conditions  for  2-D  welding  simulation  models 

The  boundary  conditions  used  in  the  model  for  the  analysis  are  shown  in  the  Figure  9.  This 
shows  the  corners  where  fillet  weld  is  done  as  heat  sources.  All  the  free  surfaces  are  taken  as  the 
convective  surfaces.  The  bottom  center  node  is  restricted  in  all  degrees  of  freedom  and  the  top 
center  node  of  the  stiffener  plate  is  restricted  in  X  direction. 

The  Finite  Element  solutions  are  performed  using  a  SMP  Fortran90  finite  element  code 
developed  in-house  for  the  2-D  models  and  ABAQUS  software  is  used  for  the  3-D  models.  The 
implementation  details  pertaining  to  those  problems  like  the  type  of  elements,  boundary  and  loading 
conditions  are  presented  to  allow  convenient  reproduction. 


2-D  Welding  Simulation 
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The  2-D  finite  element  mesh  used  in  the  heat  transfer  and  mechanical  analysis  of  the  panel  is 
illustrated  in  Figure  4.  The  model  is  made  up  of  heat  conduction,  quadratic(8  node),  quadrilateral 
elements.  The  quasi-static  mechanical  problem,  following  the  heat  transfer  analysis,  is  discretized 
into  a  generalized  plane  strain  finite  element  model.  Model  for  case  1  uses  1123  nodes  and  316 
elements.  Models  for  case  2  and  case  3  use  2973  nodes  and  3253  nodes  respectively  and  866  elements 
and  946  elements  respectively. 

To  model  the  restriction  of  the  supporting  plate,  the  downward  motion  of  the  plates  should  be 
restrained.  This  is  implemented  by  placing  two  nonlinear  spring  elements  at  the  two  outermost 
nodes  of  the  2-D  model,  to  exert  high  reaction  forces  to  resist  the  downward  motion.  The  spring 
element  stiffness  is  defined  as  a  nonlinear  function  of  the  y-displacement  as  follows, 

*  =  {“’  »>0  [N/mml  (12) 


3-D  Buckling  Analysis 


The  3-D  continuum  finite  element  models  are  developed  to  calculate  critical  buckling  loads  A  . 
2-D  generalized  plain  strain  models  are  developed  to  compute  the  stress  resulting  from  unit  thermal 
load  CL¬ 
AW  3-D  models  contain  1020  Hex  8  elements,  and  a  total  of  2184  nodes.  The  first  2-D  generalized 
plain  strain  models  contains  104  Quad  8  elements  and  423  nodes,  the  second  contains  216  Quad  8 
elements  and  767  nodes,  and  the  third  contains  216  Quad  8  elements  and  767  nodes. 

Unit  thermal  load  is  added  along  the  crossing  of  the  flanges  and  the  stiffener,  no  preload  is 
added.  The  3-D  finite  element  mesh  of  the  small  stiffener  is  illustrated  in  Figure  5. 


3,1.4  Results 


Numerical  Results 


2-D  Analysis 


After  running  the  welding  simulation  for  2-D  model,  the  distortion,  Cauchy’s  stress  and  plastic 
strain  in  longitudinal  direction  are  obtained.  Figure  10  shows  the  distribution  of  residual  stresses 
in  Z  direction  after  the  2-D  analysis  of  welding  for  case  1.  The  stresses  are  tensile  in  the  weld  zone 
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Figure  10:  Residual  Stress  obtained  from  2-D  analysis  for  case  1 


Case 

Size 

in* 

Thickness 

in. 

Heat  input 
kJ/in. 

Sres 

MPa 

i 

10x10 

1/16 

4.8 

-127 

2 

10x10 

1/8 

14.3 

-190 

3 

12x12 

1/8 

14.3 

-156 

Table  2:  Analytical  Residual  stress  using  2-D  analysis 


and  in  heat  affected  zone.  However,  away  from  the  weld  zone,  the  stress  changes  the  sign  from 
tensile  to  compressive.  Figures  11  and  12  show  the  plastic  strain  and  angular  distortion  for  case 
1.  Plastic  strain  also  changes  sign  farther  from  the  weld  zone.  Spring  elements  at  the  ends  of  the 
panel  do  not  allow  downward  distortion  thus  correctly  modelling  the  base  support  plate. 

Table  2  lists  the  values  of  maximum  compressive  stress  in  all  the  panels.  The  residual  stress  in 
the  10”  X  10”  thin  (1/16”)  panel  is  the  minimum  and  that  in  the  10”  X  10”  thick  (1/8”)  panel  is 
the  maximum.  The  table  also  lists  the  heat  input  that  was  used  to  cause  that  residual  stress.  A 
smaller  weld  is  needed  to  weld  the  thin  panel  (1/16”  thickness).  And  more  heat  input  will  result 
in  a  bum-through  of  the  specimen.  This  is  the  maximum  heat  input  that  can  be  safely  induced  in 
a  double  fillet  weld  configuration  with  a  satisfactory  fillet  weld  and  without  getting  burn-through. 
On  the  other  hand,  in  the  thick  panel,  welding  heat  is  dictated  by  the  satisfactory  fillet  weld.  Both 
the  geometries  of  the  thick  panel  use  the  same  heat  input. 
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Figure  12:  Angular  Distortion  obtained  from  2-D  analysis  for  case  1 
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Case 

Size 

in. 

Thickness 

in. 

Heat  input 
kJ/in. 

Sres 

MPa 

$cr 

MPa 

Buckling 

Predicted 

Whether 

buckled 

i 

10x10 

1/16 

4.8 

-127 

-45 

Yes 

Yes 

2 

10x10 

3/16 

14.3 

-190 

-191 

No 

No 

3 

12x12 

3/16 

14.3 

-156 

-134 

Yes 

Yes 

Table  4:  Comparison  Table  for  Analytical  and  experimental  results 


three  panels  are  welded  and  let  cool  down.  After  the  panels  are  cooled  down,  the  distortion  is 
examined.  The  10”  X  10”  X  1/16”  (case  1)  panel  shows  clear  buckling  distortion.  This  is  shown 
in  Figure  14.  However,  10”  X  10”  X  1/8”  (case  2)  and  12”  X  12”  X  1/8”  (case  3)  panels  behave 
differently  with  the  same  welding  conditions  and  hence  the  heat  input.  The  case  2  panel  shows 
angular  distortion  and  there  is  no  buckling  distortion  in  that  panel  as  shown  in  Figure  15.  On  the 
other  hand,  case  3  panel  shows  buckling  distortion  clearly.  This  buckling  can  be  seen  in  Figure  16. 


Figure  14:  Buckling  Distortion  in  case  1 

A  table  to  compare  2-D  and  3-D  analytical  results  with  experimental  results  is  drawn.  Table 
4  shows  the  successful  validation  of  analytical  results.  The  experimental  results  are  in  agreement 
with  the  corresponding  analytical  results.  In  the  cases  where  buckling  is  predicted  analytically,  we 
can  see  the  occurrence  of  buckling  in  the  experimental  results. 


3.2  Inverse  method  and  database  architecture 


Sensitivity  analysis  has  been  developed  for  the  eigenvalue  buckling  problem.  The  derivative  of 
a  panel’s  critical  buckling  strength  with  respect  to  structural  parameters  such  as  panel  thickness 
and  size  is  computed  in  numerically  consistent  and  efficient  manner.  This  section  presents  the 
mathematical  formulations  and  verification  performed  on  a  small  example. 
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Figure  16:  Buckling  Distortion  In  Case  3 


3.2.1  Mathematic  formulation 


Buckling  equations 


Buckling  analysis  can  be  expressed  as  a  generalized  eigenvalue  problem  as  follows  [25,  26]: 


[K  -  A»K =  0  (13) 

where  K  is  linear  stiffness  matrix,  which  is  symmetric  positive  definite.  Ks  is  stress  stiffness 
matrix  with  unit  load,  which  is  symmetric,  but  not  necessary  positive  definite,  A  is  eigenvalue, 
which  represents  the  critical  buckling  load,  and  <f>  is  eigenvector,  which  corresponds  to  buckling 
mode. 

The  detailed  formulation  of  K  and  Ks  are  given  below: 

K  =  jBTBBdV  e  =  BU  U  =  Nu  (14) 


where  N  is  the  matrix  of  shape  function,  E  is  the  matrix  of  Young’s  modula,  and  B  is  the 
derivative  of  matrix  N,  B  =  <9N.  p  is  surface  pressure  term,  and  b  is  body  force. 


Sensitivity  equations 


Sensitivity  analysis  is  one  of  the  important  aspects  in  engineering  design  theory  and  application. 
The  expressions  of  eigenvalue  sensitivity  and  eigenvector  sensitivity  are  computed  by  deriving  the 
first  order  sensitivities  for  symmetric  positive  definite  eigenvalue  systems  (13).  The  equations  are 
given  by  D.  A.  Tortorelli  and  P.  Michaleris  [27]: 


Ddi 


pK  N  DEM 
Ddi  Aj  Ddi . 


+i 


(16) 
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where  cfi®  is  the  normalized  eigenvectors  of  4>},  and  d  is  design  parameter  vector. 


The  derivative  of  eigenvectors  are  given  by  the  equation: 


[K  -  AjK3] 


Ddi 


DK  DKS' 

Ddi  Xj  Ddi . 


+J- 


D>± 

Dxi 


(17) 


Noticing  [K  -  AjKs]  is  singular  now,  first  assume  the  s-th  component  4>js  =  1,  thus  =  0, 
remove  the  corresponding  row  and  column  in  [K  -  Aj(x)Ks],  then  comes  a  equation  from’ which 

other  components  of  can  be  solved. 


3.2.2  Numerical  implementation  and  verification 

The  buckling  and  sensitivity  equations  has  been  integrated  into  FORTRAN  code,  and  used  in 
the  following  numerical  examples.  Results  are  presented  here  for  two  simple  examples.  The 
analytical  solution  along  with  results  obtained  using  the  commercial  code  ABAQUS  are  presented 
for  verification. 


Model  information 


Figure  17:  cantilever  beam,  1st  buckling  mode 
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Figure  18:  load  information 


The  first  example  is  a  cantilever  beam,  all  three  degrees  of  freedom  for  the  nodes  at  one  end  and 
one  vertical  degree  of  freedom  for  nodes  at  the  other  end  are  constrained.  A  unit,  surface  pressure 
is  applied  at  the  latter  surface  as  load.  Please,  refer  to  Figure  17  for  the  model  information  and 
Figure  18  for  the  load  information. 

The  second  example  is  a  simple  supported  beam,  where  at  the  fixed  end  all  three  degrees  of 
freedom  are  constrained  and  a  unit  surface  pressure  is  applied  at  the  other  end.  See  Figure  19  for 
the  model  information.  The  load  information  is  the  same  as  the  first  one. 

The  length,  width  and  height  of  both  beams  are  10,  1,  1  respectively,  with  Young’s  modula 
E  =  2.0E5  and  Poisson’s  ratio  7  =  0. 


Results  and  comparison 


Table  3.2.2  presents  a  comparison  of  the  eigenvalues  for  the  two  example  problems.  Analytical 
results,  computational  results  using  the  commercial  code  ABAQUS  and  the  in-house  FORTRAN 
code  are  listed.  Table  3.2.2  presents  a  comparison  of  the  sensitivities  of  eigenvalues  with  respect 
to  the  beam  length  for  the  two  example  problems.  Analytical  results,  computational  results  using 
the  direct  differentiation  and  finite  difference  of  the  in-house  FORTRAN  code  are  listed. 

For  the  cantilever  beam,  the  approximate  theoretical  critical  buckling  load  can  be  expressed  as: 
Per  =  =  411.23,  where  <f)  =  4>csin( |g)  is  used  as  the  approximate  buckling  mode.  The  result 

from  ABAQUS  is  Pcr  —  —889.21,  and  from  the  FORTRAN  code  is  P„  =  460.04.  The  buckling 
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Figure  19:  simple  supported  beam,  1st  buckling  mode 


modes  are  same  from  ABAQUS  and  the  FORTRAN  code,  as  in  Figure  17. 

The  sensitivity  of  the  eigenvalue  with  respect  to  the  length  of  the  beam  as  calculated  from 
FORTRAN  code  using  Equation  (16)  is  A P  =  -80.3.  To  verify  this  result,  the  sensitivity  is  also 
computed  by  the  finite  difference  method.  The  length  of  beam  is  changed  to  10.1,  with  other 
properties  remaining  the  same,  the  critical  buckling  load  changes  to  =  452.04  which  results  to 
a  finite  difference  sensitivity  of  80.0. 

For  the  simply  supported  beam,  the  approximate  theoretical  critical  buckling  load  is  P„  = 
3^1  =  1644.93,  where  <f>  =  <f>csin( is  used  as  approximate  buckling  mode.  The  result  from 
ABAQUS  is  Pcr  =  1633.7,  and  from  FORTRAN  code  is  Per  =  1688.81.  The  result  of  buckling 
modes  are  also  same  from  ABAQUS  and  FORTRAN  code,  as  in  Figure  19. 

For  the  finite  difference  sensitivity,  the  critical  buckling  load  for  the  beam  with  length  10.1  is 
Per  —  1661.58  which  results  to  the  finite  difference  sensitivity  of  A P  —  -272.3.  The  numerical 
sensitivity  result  calculated  from  FORTRAN  code  is  A P  =  —272.8. 


model 

Analytical  method 

ABAQUS 

FORTRAN  code 

i 

411.23 

-889.21 

460.04 

2 

1644.93 

1633.7 

1688.81 

Table  5:  Eigenvalues  calculated  from  analytical  method,  ABAQUS  and  FORTRAN  code 
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model 

Analytical  method 

Finite  difference  method 

FORTRAN  code 

i 

-82.2 

-80.0 

-80.3 

2 

-329.0 

-272.3 

-272.8 

Table  6:  Sensitivities  calculated  from  analytical  method,  finite  difference  and  FORTRAN  code 


3.3  Generation  of  a  database 


Work  in  this  task  will  commence  in  FY  02. 


3.4  Integration  with  computer  aided  production  engineering  software 

Work  in  this  task  will  commence  in  FY  02. 


3.5  Modifications  and  refinement  using  production  tests 


Work  in  this  task  will  commence  in  FY  02. 
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